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ABSTRACT 

This  document  describes  a preliminary  version  of  a long  period  array 
processing  package  designed  around  the  FKCOMB  algorithm  for  use  on  the 
ILLIAC  IV  computer.  FKCOMB  is  a general-purpose  array-processing  program 
that  uses  frequency-wavenumber  analysis  to  produce  a bulletin  which  lists 
signal  detections  and  various  statistics  for  each  detection.  Two  data 
editing  and  reformatting  modules  prepare  the  seismic  data  for  FKCOMB  and 
can  be  modified  for  use  with  other  seismic  algorithms.  Preliminary  ref or- 
matting  of  the  seismic  data  is  performed  "by  DEMI.  The  data  is  edited  and 
fast  fourier  transformed  by  DEM2. 

The  input  parameters  required  for  operating  these  programs  and  their 
subroutines  are  described  in  this  document. 
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INTRODUCTION 


The  preliminary  version  of  the  long  period  array  processing  package  is 
designed  to  demonstrate  the  feasibility  of  using  the  ILLIAC  IV  computer  for 
seismic  processing.  It  consists  of  three  modules.  Data  editing  module  one 
(DEMI)  reformats  the  long  period  data,  Data  editing  module  two  (DEM2)  edits 
and  fast  Fourier  transforms  (FFT)  the  data  and  FKCOMB  performs  the  FKCOMB 
algorithm  on  the  data.  Each  of  these  modules  is  structured  to  allow  for 
expansion  to  include  additional  data  formats. 

Figure  1 is  an  overall  flow  of  data  from  SDAC  to  the  ILLIAC  site. 
Currently  the  data  is  received  at  ILLIAC  over  the  ARPANET  via  file  transfer 
protocol  (FTP)  in  the  SDAC  low-rate  tape  format.  A tape  drive  is  available 
at  the  ILLIAC  site  and  could  be  used  for  data  transfer.  Output  from  the 
program  is  transmitted  to  SDAC  via  the  ARPANET. 

The  channel  component  data  is  retrieved  by  DEMI  and  stored  by  array 
identifier  in  time  order.  Minor  timing  errors  are  corrected.  DEM2  edits 
the  data  for  timing  and  data  errors,  transforms  the  data  and  stores  the 
output  in  a disk  file  for  input  to  the  FKCOMB  algorithm. 
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DATA  EDITING  MODULE  1 (DEMI) 


PURPOSE 


Seismic  data  as  received  from  remote  seismic  arrays  (LASA,  NORSAR,  ALPA) 
is  not  conveniently  formatted  for  efficient  retrieval  by  detection  and  analy- 
sis programs.  As  received,  records  from  different  arrays  are  randomly 
intermixed.  The  format  and  length  of  each  record  is  dependent  upon  its 
site  or  origin.  The  density  of  data  pertinent  to  most  programs  is  under 
20%.  Over  90%  of  the  runtime  of  the  current  version  of  FKCOMB  is  spent 
getting  the  data  in  a useful  form.  This  problem  is  heightened  by  the 
parallel  architecture  of  ILL1AC.  In  its  raw  form  seismic  data  would  be  very 
difficult  to  manipulate  in  a parallel  fashion.  DEMI  is  designed  to  read  in 
raw  data  and  organize  it  in  a convenient  manner  for  processing  by  such 
algorithms  as  FKCOMB.  The  structure  of  the  data  makes  it  difficult  to 
maximize  the  processing  power  available  on  ILLIAC  during  such  a reorganiza- 
tion. This  is  not  important,  as  the  process  is  inherently  I/O  bound. 
Utilization  of  the  bandwidth  between  disk  and  core  (,5*10**9  bits/sec.)  and 
the  large  amount  of  core  available  (128K  64-bit  words)  makes  it  possible  to 
do  this  job  on  ILLIAC  quickly  and  makes  the  data  available  for  storage  on 
the  UNICON  laser  memory  for  convenient  access  for  any  processing  desired. 

The  program  is  flexible  and  could  handle  data  from  all  existing  arrays  and 
from  any  similar  long  period  or  short  period  arrays  available  in  the  future. 


FUNCTIONAL  AND  THEORETICAL  DISCUSSION 


The  input  and  output  files  used  by  DEMI  may  be  thought  of  as  a series 
of  16-bit  bytes.  It  is  the  function  of  DEMI  to  retrieve  the  data  bytes 
from  the  input  file  and  organize  them  into  the  appropriate  location  in  the 
output  file.  Other  than  some  minor  error  correction,  the  order  of  the  input 
bytes  and  the  bytes  themselves  remain  unchanged.  Since  the  format  of  the 
input  records  is  known,  the  task  reduces  to  determining  which  byte  of  the 
input  file  is  to  be  used  to  fill  each  byte  of  each  output  file,  and  having 
determined  this  to  access  the  bytes  in  such  a way  to  minimize  time  lost  due 


-3- 


to  memory  and  disk  accesses.  In  practice,  the  method  is  turned  around 
somewhat.  Since  the  order  of  the  input  bytes  is  unchanged,  it  is  effective 
to  take  bytes  from  the  input  file  sequentially  and  insert  each  one  into  the 
next  location  in  the  appropriate  output  file.  In  order  to  conserve  disk 
access  time,  large  core  buffers  are  used  for  input  and  output.  ILLIAC  core 
contains  approximately  8x10^  bits.  Using  buffers  of  10^  bits,  500  disk 

g 

accesses  are  necessary  to  input  5x10  bits  (24  hours  of  data)  and  another 
500  disk  accesses  to  output  the  data.  Assuming  the  worst  case  of  40  milli- 
second access  time  this  means  that  approximately  40  seconds  of  I/O  time  will 
be  required  to  process  24  hours  of  data  using  a simple  single  buffering 
scheme  for  1/0,  if  all  of  core  is  used  for  buffering.  In  actuality,  not  all 
of  core  is  available  for  1/0  buffers  and  the  time  taken  may  be  up  to  80 
seconds.  Experimentation  with  different  disk  mappings  to  reduce  the  average 
access  time  may  significantly  reduce  the  time.  The  second  major  source  of 
time  is  memory  access.  ILLIAC  does  not  allow  memory  to  memory  transfers, 
so  all  data  most  go  from  memory  to  a register  to  memory.  ILLIAC  memory 
access  time,  though  variable,  due  to  overlap,  is  approximately  600  nano- 
seconds. If  a transfer  were  done  for  each  16-bit  byte,  requiring  two  memory 
accesses,  2 accesses  x (6x10  ) seconds/access  x 10  bytes  = 12  seconds  of 

memory  access  time  would  be  required  to  transfer  all  input  bytes.  As  coded, 
the  scratch  pad  memory  in  the  ILLIAC  Control  Unit  (CU)  is  used  to  somewhat 
reduce  this  time. 

The  third  major  source  of  time  is  overhead  due  to  the  calculation  of 
the  address  in  each  buffer  before  a transfer  and  the  shifting  necessary  to 
coordinate  the  16-bit  byte  size  of  the  data  with  the  64-bit  ILLIAC  word 
size.  The  amount  of  time  used  in  this  process  is  estimated  to  be  equal  to 
or  slightly  greater  than  the  time  taken  in  memory  accesses. 

The  time  required  by  DEMI  to  reorganize  twenty  four  hours  of  raw  long 
period  seismic  data  is  under  two  minutes  of  1/0  time  and  under  one  minute 
of  processor  time.  Allowing  for  unforeseen  overhead  and  future  code 
optimization,  three  to  five  minutes  of  total  running  time  is  a reasonable 
estimate  for  this  system  of  data  acquisition.  More  complicated  and  possibly 
more  efficient  algorithms  suggest  themselves  but  are  deemed  unnecessary  due 
to  the  small  savings  in  time  possible. 

-4- 


r-'-U bVi vK’ff/rir- , 


mmm 


PROG  RAM  DESCRIPTION 


[ 

Seismic  data,  as  received  over  remote  lines  is  divided  into  records. 
Each  record  consists  of  a series  of  16-bit  bytes.  While  the  Length  and 
exact  format  of  each  record  depends  upon  the  array  at  which  it  was  created, 
the  general  format  of  each  is: 


Record  Format 
HEADER  ID  - 

POSSIBLE  TIMING  WORD  - 

VARIABLE  LENGTH  GAP  - 
MULTIPLE  TIME  SCANS  - 

VARIABLE  LENGTH  GAP  - 
Time  SCAN 

VARIABLE  LENGTH  GAP- 
POSSIBLE  TIME  WORD  - 


MULTTPLE  DATA  CHANNELS  - 
VARIABLE  LENGTH  GAP  - 


Identifies  source  array. 

Present  if  one  time  word  per  record 
rather  than  one  time  word  per  scan  (for 
definition  of  SCAN,  see  below).  Format 
of  time  word  dependent  upon  array. 

Length  determined  by  array. 

Number  dependent  upon  array.  For  format 
of  SCAN,  see  below. 

Length  determined  by  array. 

Length  determined  by  array. 

Present  if  time  word  per  SCAN  rather 
than  time  word  per  record.  Format  of 
time  word  dependent  upon  array. 

Number  and  format  dependent  upon  array. 
Length  determined  by  array. 


In  this  general  structure  there  are  10  variables  which  exactly  determine 
the  format  of  a record.  They  are: 

1)  HEADER  ID 

2)  TIME  WORD  PER  SCAN  OR  TIME  WORD  PER  RECORD 

3)  FORMAT  OF  TIME  WORD 

4)  LENGTH  OF  GAP  BEFORE  FIRST  TIME  SCAN 

5)  LENGTH  OF  GAP  AT  BEGINNING  OF  EACH  SCAN 

6)  NUMBER  OF  DATA  CHANNELS  PER  SCAN 
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7)  FORMAT  OF  DATA  CHANNELS 

8)  LENGTH  OF  GAP  AT  END  OF  EACH  SCAN 

9)  LENGTH  OF  GAP  AFTER  LAST  SCAN 

10)  NUMBER  OF  SCANS  PER  RECORD. 

DEMI  determines  the  header  ID  from  the  first  byte  of  a record  and  from 
this  determines  the  source  array.  All  the  other  values  are  then  found  in  a 
table  constructed  to  describe  each  array.  New  formats  can  be  handled  by 
adding  to  or  modifying  this  table.  The  table  is  currently  constructed  at 
compile  time  via  a data  statement  and  is  changed  by  recompiling  one  block 
data  subroutine. 

A pointer  is  maintained  indicating  which  byte  in  the  input  stream  is  to 
be  accessed  next.  Using  this  pointer,  it  is  determined  whether  the  data 
byte  is  on  disk,  in  core,  or  in  the  CU  scratch  pad  memory.  The  format  of 
the  input  is  accommodated  by  manipulating  the  pointer  to  space  over  gaps. 

In  a series  of  nested  loops,  each  record,  then  each  time  scan,  then  each 
data  channel  is  handled.  Output  goes  to  one  of  six  output  buffers,  then  to 
one  of  six  output  files.  Since  only  three  seismic  arrays  are  currently 
implemented,  only  three  of  these  files  are  presently  used.  The  others  are 
available  for  expansion. 

Timing  errors  in  the  input  are  of  two  types,  gaps  and  reversals.  Gaps 
of  two  or  fewer  are  corrected  by  duplicating  the  data  from  the  last  good 
time  step.  This  data  is  always  available  in  the  output  buffer,  though  the 
duplication  is  complicated  somewhat  by  the  two  levels  of  buffering  used. 
Larger  gaps  are  indicated  by  a diagnostic  message  and  no  correction  is 
attempted.  Time  reversals  are  handled  by  ignoring  any  redundant  data  and 
a diagnostic  message. 

Processing  continues  until  an  unidentifiable  input  record  is  found. 

If  the  header  is  non-zero,  a diagnostic  is  printed  indicating  the  possibility 
of  a data  error.  Normal  end  of  data  is  indicated  by  a header  of  all  zeroes. 
In  either  case,  all  output  buffers  are  emptied  to  disk  and  processing  ter- 
minates . 


DATA  AREAS  AND  SYMBOL  DEFINITIONS 
ADB  - CU  INTEGER 

ADBBUF(8)  - CU  INTEGER 

AD BOUT (6)  - CU  INTEGER 
ADDRS  - CU  INTEGER 
ARRAY  - CU  INTEGER 
BCT  - CU  INTEGER 

BYTCNT(6)  - CU  INTEGER 
BYTS  - CU  INTEGER 


CNTRL  (*,6)  _ PE  INTEGER 

DEBUG  - CU  INTEGER 

END  ADB  - CU  INTEGER 

INBUF  (* ,128)  - PE  INTEGER 
TNBYT  - CU  INTEGER 
INPTB  - CU  INTEGER 
INPTW  - CU  INTEGER 
IT  - CU  INTEGER 

LADB  - CU  LOGICAL 

LADBBU (8)  - CU  LOGICAL 

LADBOU  - CU  LOGICAL 


- holds  word  from  ADB  buffer  before 
being  written  in  core. 

- CU  scratch  pad  memory  (ADB)  input 
buffer. 

- one  word  ADB  buffer  for  each  array. 

- address  in  OUTBUF  of  IT. 

- Array  currently  being  processed. 

- number  of  bytes  in  last  partially 
filled  word. 

- byte  count  for  each  output  buffer. 

- Number  of  bytes  from  last  time  step 
that  are  in  core.  Used  in  filling 
time  gaps. 

- An  array  initialized  at  compile  time 
giving  the  format  of  data  from  each 
seismic  array. 

- Controls  the  output  of  debug  print- 
out. Zero  for  typical  run. 

- address  of  byte  being  held  at  end  of 
ADB  input  buffer. 

- Input  buffer. 

- holds  input  byte  after  call  to  GETBYT. 

- Pointer  to  current  byte  in  input  buffer. 

- Pointer  to  current  word  in  input  buffer. 

- first  word  of  data  to  use  when  filling 
in  time  gaps. 

- equivalenced  to  ADB  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ADBBUF  to  facilitate 
shifting  and  masking. 

- equivalencea  to  ADBOUT  to  facilitate 
shifting  and  masking. 


LADBWR  - CU  LOGICAL 
LADDRS  - CU  LOGICAL 
LARRAY  - CU  LOGICAL 
LBCT  - CU  LOGICAL 
LBYTCN(6)  - CU  LOGICAL 
LBYTS  - CU  LOGICAL 
LDEBUG  - CU  LOGICAL 
LENDAD  - CU  LOGICAL 
LINBYT  - CU  LOGICAL 
LINPTB  - CU  LOGICAL 
LINPTW  - CU  LOGICAL 
LIT  - CU  LOGICAL 
LORGCO  - CU  LOGICAL 
LOUBYT  - CU  LOGICAL 
LOUPTW  - CU  LOGICAL 
LPAGE  - CU  LOGICAL 
LT1  - CU  LOGICAL 


- equivalenced  to  ADBWRD  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ADDRS  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ARRAY  to  facilitate 
shifting  and  masking. 

- equivalenced  to  BCT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  BYTCNT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  BYTS  to  facilitate 
shifting  and  masking. 

- equivalenced  to  DEBUG  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ENDADB  to  facilitate 
shifting  and  masking. 

- equivalenced  to  INBYT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  INPTB  to  facilitate 
shifting  and  masking. 

- equivalenced  to  INPTW  to  facilitate 
shifting  and  masking. 

- equivalenced  to  IT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ORGCOR  to  facilitate 
shifting  and  masking. 

- equivalenced  to  OUBYT  co  facilitate 
shifting  and  masking. 

- equivalenced  to  OUPTW  to  facilitate 
shifting  and  masking. 

- equivalenced  to  PAGE  to  facilitate 
shifting  and  masking. 

- equivalenced  to  TI  to  facilitate 
shifting  and  masking. 


DATA  AREAS  AND  SYMBOL  DEFINITIONS 
ADB  - CU  INTEGER 

ADBBII  > - CU  INTEGER 

AD BOUT (6)  - CU  INTEGER 
ADDRS  - CU  INTEGER 
ARRAY  - CU  INTEGER 
BCT  - CU  INTEGER 

BYTCNT (6)  - CU  INTEGER 
BYTS  - CU  INTEGER 


CNTRL  (*,6)  _ PE  INTEGER 

DEBUG  - CU  INTEGER 
END  ADB  - CU  INTEGER 
INEUF  (* ,128)  - PE  INTEGER 

£ 

INBYT  - CU  INTEGER 
INPTB  - CU  INTEGER 
INPTW  - CU  INTEGER 
IT  - CU  INTEGER 

LADB  - CU  LOGICAL 

LADBBU (8)  - CU  LOGICAL 

LADBOU  - CU  LOGICAL 

1 

I 
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- holds  word  from  ADB  buffer  before 
being  written  in  core. 

- CU  scratch  pad  memory  (ADB)  input 
buffer. 

- one  word  ADB  buffer  for  each  array. 

- address  in  OUTBUF  of  IT. 

- Array  currently  being  processed. 

- number  of  bytes  in  last  partially 
filled  word. 

- byte  count  for  each  output  buffer. 

- Number  of  bytes  from  last  time  step 
that  are  in  core.  Used  in  filling 
time  gaps. 

- An  array  initialized  at  compile  time 
giving  the  format  of  data  from  each 
seismic  array. 

- Controls  the  output  of  debug  print- 
out. Zero  for  typical  run. 

• address  of  byte  being  held  at  end  of 
ADB  input  buffer. 

- Input  buffer. 

- holds  input  byte  after  call  to  GETBYT. 

- Pointer  to  current  byte  in  input  buffer. 

- Pointer  to  current  word  in  input  buffer. 

- first  word  of  data  to  use  when  filling 
in  time  gaps . 

- equivalenced  to  ADB  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ADBBUF  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ADBOUT  to  facilitate 
shifting  and  masking. 
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LADBWR  - CU  LOGICAL 
LADDRS  - CU  LOGICAL 
LARRAY  - CU  LOGICAL 
LBCT  - CU  LOGICAL 
LBYTCN(6)  - CU  LOGICAL 
LBYTS  - CU  LOGICAL 
LDEBUG  - CU  LOGICAL 
LENDAD  - CU  LOGICAL 
LINBYT  - CU  LOGICAL 
LINPTB  - CU  LOGICAL 
LINPTW  - CU  LOGICAL 
LIT  - CU  LOGICAL 
TORGCO  - CU  LOGICAL 
LOUBYT  - CU  LOGICAL 
LOUPTW  - CU  LOGICAL 
LPAGE  - CU  LOGICAL 
LT1  - CU  LOGICAL 


- equivalenced  to  ADBWRD  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ADDRS  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ARRAY  to  facilitate 
shifting  and  masking. 

- equivalenced  to  BCT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  BYTCNT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  BYTS  to  facilitate 
shifting  and  masking. 

- equivalenced  to  DEBUG  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ENDADB  to  facilitate 
shifting  and  masking. 

- equivalenced  to  INBYT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  INPTB  to  facilitate 
shifting  and  masking. 

- equivalenced  to  INPTW  to  facilitate 
shifting  and  masking. 

- equivalenced  to  IT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  ORGCOR  to  facilitate 
shifting  and  masking. 

- equivalenced  to  OUBYT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  OUPTW  to  facilitate 
shifting  and  masking. 

- equivalenced  to  PAGE  to  facilitate 
shifting  and  masking. 

- equivalenced  to  TI  to  facilitate 
shifting  and  masking. 
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LT2  - CU  LOGICAL 

equivalenced  to  T2  to  facilitate 
shifting  and  masking. 

LT3  - CU  LOGICAL 

“* 

equivalenced  to  T3  to  facilitate 
shifting  and  masking. 

LT4  - CU  LOGICAL 

— 

equivalenced  to  T4  to  facilitate 
shifting  and  masking. 

LT5  - CU  LOGICAL 

equivalenced  to  T5  to  facilitate 
shifting  and  masking. 

LT6  - CU  LOGICAL 

equivalenced  to  T6  to  facilitate 
shifting  and  masking. 

LWORD  - CU  LOGICAL 

equivalenced  to  WORD  to  facilitate 
shifting  and  masking. 

LWORDS  - CU  LOGICAL 

equivalenced  to  WORDS  to  facilitate 
shifti.  g and  masking. 

LPRTIA  - CU  LOGICAL 

equivalenced  to  PARTIAL  to  facilitate 
shifting  and  masking. 

LSAVAD  - CU  LOGICAL 

“ 

equivalenced  to  SAVADB  to  facilitate 
shifting  and  masking. 

OLDTIM(*)  - PE  INTEGER 

Time  in  deciseconds  from  the  beginning 
of  the  year  for  current  time  step. 

ORGCOR  - CU  INTEGER 

identifies  address  within  file  of  the 
first  word  in  the  core  input  buffer. 

OTIMEA(6)  - PE  INTEGER 

Saves  OLDTIM  when  switching  from  one 
array  to  another. 

OUBYT  - CU  INTEGER 

- 

passes  output  byte  to  subroutine  PUTBYT 

OUPAGE(6)  - PE  INTEGER 

- 

page  of  current  output  file 

OUPTWA(6)  - PE  INTEGER 

Pointer  to  current  word  in  output  buffe: 
for  each  array. 

OUT(*, 64,6)  - PE  INTEGER 

Output  buffer.  Includes  room  for  6 
seismic  arrays. 

OUTPTW  - CU  INTEGER 

pointer  to  current  word  in  output 
buffer. 

PAGE  - CU  INTEGER 

- 

page  to  read  from  next. 
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PINTM*)  - PE  INTEGER 

PRTIAL  - CU  INTEGER 

PaVABD  - CU  INTEGER 

SAVBCT  - PE  INTEGER 

SAVPTW  - PE  INTEGER 
SCANS  - PE  INTFGER 

TD1E(*)  - PE  INTEGER 

TSTEPS (6)  - PE  INTEGER 

T1  - CU  INTEGER 

T2  - CU  INTEGER 

T3  - CU  INTEGER 

T4  - CU  INTEGER 

T5  - CU  INTEGER 

T6  - CU  INTEGER 

WORD  - CU  INTEGER 

WORDS  - CU  INTEGER 


- Holds  intermediate  results  throughout 
the  program. 

- Number  of  bits  of  IT  to  use  when 
filling  in  time  gaps. 

- Saves  ADBOUT( ARRAY)  while  filling  in 
time  gaps. 

- Saves  BYTCNT  (ARRAY)  while  filling  in 
time  gaps. 

- Saves  OUPTW  while  filling  in  time  gaps. 

- Number  of  time  scans  that  have  been 
read  from  the  current  record. 

- Time  in  deciseconds  from  the  beginning 
of  the  year  for  current  time  step. 

- Number  of  time  steps  that  have  been 
retrieved  for  each  array. 

- holds  intermediate  results  throughout 
program. 

- holds  intermediate  results  throughout 
program. 

- holds  intermediate  results  throughout 
program. 

- holds  intermediate  results  throughout 
program. 

- holds  intermediate  results  throughout 
program. 

- holds  intermediate  results  throughout 
program. 

- used  to  build  output  word  while  filling 
time  gaps. 

- number  of  complete  words  that  are  in 
core  from  last  time  window.  Used  in 
filling  time  gaps. 


LINKAGE 


J 

DEMI  calls  subroutine  PUTBYT , GETBYT,  CNVTIM,  and  RDPRM  which  are 
described  in  the  section  on  subroutines . 

DISK  AREAS 

INPUT  - binary  input  area.  Format  described  in  DEMI  program 
description. 

OUPUT1,  0UPUT2 , 0UPUT3,  0UPUT4,  0UPUT5,  OUPUT6  - binary  output  areas. 

Each  contains  output  from  one  seismic  array.  Currently  the  input  consists 
of  data  from  only  three  arrays,  so  only  the  first  three  areas  contain  data 
other  than  the  header.  The  format  of  each  area  is: 

HEADER  - beginning  of  each  area. 

WORDl:  Array  ID. 

W0RD2 : number  of  time  steps  in  this  disk  area. 

WORD3-1024 : zero 

TIME  STEPS  - first  time  step  starts  at  word  1025.  Number  of  time 

5 

steps  given  in  word  2 of  the  HEADER. 

WORDl:  time  in  deciseconds  from  beginning  of  year  for  this 

time  step. 

WORD2-N : data  from  each  of  N— 1 channels. 

DISPA  - Print  output.  Consists  of  a header  identifying  the  run, 
bulletins  indicating  time  gaps  and  reversals,  and  a trailer  indicating 
normal  termination  of  the  run. 


'm 


DATA  EDITING  MODULE  2 (DEM2) 


PURPOSE 

The  FKCOMB  algorithm  requires  that  the  input  data  be  rearranged,  edited, 
FFT  ed,  and  rearranged  again  before  detection  can  begin.  For  ease  of  design 
and  because  of  core  limitations,  this  part  of  the  processing  has  been  made 
a separate  module.  DEM2  accepts  the  output  of  DEMI  and  creates  an  output 
file  which  is  error  free  and  ready  for  FKCOMB  detection.  It  handles  any 
overlap  requested  and  prints  a bulletin  whenever  time  gaps  or  data  errors 
are  detected.  The  data  is  deglitched  and  dead  or  noisy  channels  are 
removed. 

FUNCTIONAL  AND  THEORETICAL  DISCUSSION 

Long  period  seismic  data  transmitted  to  SDAC  is  multiplexed  by  time. 

This  organization  o.s  unaffected  by  DEMI.  Prior  to  FKCOMB  analysis,  the 
following  steps  must  be  performed  on  the  data: 

• demultiplex  by  channel  and  segment  into  time  windows 

• overlap  of  time  windows  according  to  user  specification 

• convert  to  ILLIAC  floating  point  format 

• correct  large  time  gaps  by  shifting  of  time  windows 

• deglitch 

• remove  dead  and  noisy  channels 

• fast  fourier  transform 

• demultiplex  by  frequency. 

Both  steps  involving  demultiplexing  are  done  a byte  at  a time.  All 
other  steps  make  full  use  of  the  parallel  structure  of  ILLIAC  by  moving 
or  computing  64  pieces  of  data  at  once. 

The  correction  of  time  gaps  and  th?  overlapping  of  time  windows 
involves  moving  data  across  processing  elements  (PE).  Due  to  the  organi- 
zation of  data  into  time  windows  whose  size  is  a power  of  two  equal  to  or 
greater  than  64,  tnis  is  an  efficient  use  of  the  ILLIAC  route  instruction. 
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Spikes  are  smoothed  via  the  following  formula.  Let  An-1,  A and 
V-!  be  three  consecutive  data  elements  and  G be  the  user  supplied  "glitch" 
factor”.  If  Ka.-A^MgCA^  -A^)  | than  An  - (A^  + A^/2. 

Dead  and  noisy  channels  are  detected  via  the  mean  square.  If  C is  the 
mean  square  of  the  nth  channel,  MS  is  the  mean  square  of  all  the  channels 
and  V is  the  user  supplied  "variance  factor",  a check  is  made  to  see  if 
1/V  < cn^  < v*  not«  the  channel  is  assumed  to  be  dead  or  noisy. 

PROGRAM  DESCRIPTION 

Raw  seismic  data  and  data  as  output  by  DEMI  are  arranged  by  time  step 
in  the  following  manner: 

T(1)[CH(1),CH(2),...,CH(N)3,  T(2) [CH(1) ,CH(2) CH(N)] 

T(M)[CH(1),CH(2),...,CH(N)] 

where  "N"  is  the  number  of  channels  of  the  array  being  processed  and  "M" 
is  the  number  of  time  steps  contained  in  the  data  set  being  processed.  As 
output  by  DEMI,  the  data  is  packed  four  16-bit  bytes  per  ILLIAC  word.  The 
first  step  of  DEM2  is  to  read  in  a time  window  of  data  and  multiplex  it  by 
channel  so  that  it  is  arranged  as  follows: 

CH(1) [T(l) ,T(2) , . . . ,T(W) ] CH (N) [T(l) ,T(2) , . . ,T(W) ] 

where  "N"  is  again  the  number  of  channels  of  the  array  being  processed  and 
"W"  is  the  time  window  length,  a power  of  two  between  64  and  512  specified 
by  the  user  at  run  time.  During  this  process,  the  data  is  unpacked  to  one 
16-bit  byte  per  ILLIAC  word.  A double  buffering  scheme  is  used  so  that  the 
last  time  window  craeated  is  always  accessible  so  that  any  overlap  specified 
by  the  user  can  be  handled  and  any  time  gaps  recovered  from,  is  possible. 

The  data  is  then  deglitched  and  bad  channels  are  noted.  The  data  is 

then  FFT'ed.  The  FFT  routine  is  one  written  at  the  University  of  Illinois. 

It  is  necessary  to  convert  the  data  to  32-bit  floating  point  format  before 
the  routine  is  called.  The  FFT  is  performed  in  place  and  the  output  is  a 
complex  number  stored  in  32-bit  form  the  inner  and  outer  parts  of  the  64-bit 

word.  The  data  is  left  in  this  form  for  use  by  FKCOMB.  The  last  step  is  to 
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select  the  frequencies  requested,  rearrange  the  data  so  that  each  PE  memory 
contains  one  complete  multi-channel  time  window  multiplexed  by  frequency. 
640  words  in  each  PE  are  reserved  for  the  output,  so  the  number  of  channels 
multiplied  by  the  number  of  frequencies  is  restricted  to  less  than  640. 
Since  in  normal  operation  fewer  than  20  frequencies  are  requested,  this 
poses  no  serious  restriction. 


DEM2  VARIABLES 

ABUFF2( 70400)  - PE  INTEGER 
ADBBUF(8)  - CU  INTEGER 
ADBWRD  - CU  INTEGER 
ALLMSQ(*)  - PE  REAL 

ARRAY  - CU  INTEGER 


BF3PE  - CU  INTEGER 

BUFF2(*,500,2)  - PE  REAL 

BUFF3(8,640)  - PE  REAL 
BYTE  - CU  INTEGER 

CH  - CU  INTEGER 
CHG00D(80)  - PE  INTEGER 

CHMSQ( 80)  - PE  REAL 
CNTRL(*,6)  - PE  INTEGER 


C0MP(*)  - PE  INTEGER 


- equivalenced  to  BUFF2. 

- ADB  input  buffer. 

- current  word  in  ADBBUF. 

- the  metn  square  of  all  channels  for 
current  time  window. 

- indicates  seismic  array  being  pro- 
cessed. 1 = LASA 

2 = ALPA 

3 = NORSAR. 

- PE  to  be  filled  next  with  a time 
window. 

- intermediate  data  buffer  in  which 
time  windows  are  built. 

- output  buffer. 

- byte  within  current  ADB  word  most 
recently  acquired. 

- channel  currently  being  processed. 

- indicates  which  channels  passed 
variance  test. 

- mean  square  for  each  channel. 

- via  assembled  in  data,  give  informa- 
tion on  each  array  such  as  how  many 
sensors. 

- component  of  motion  to  be  processed. 

0 - vertical 

1 - horizontal. 


COREPT  - CU  INTEGER 
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DEBUG  - CU  INTEGER 


DIFFR  - CU  INTEGER 


DIFFW  - CU  INTEGER 

F - CU  INTEGER 
FINSCN(*)  - PE  INTEGER 
GAP  - CU  INTEGER 

GLCHFT(*)  - PE  REAL 

HI  FREQ  - PE  INTEGER 
I BUFFI (4096)  - PE  INTEGER 
IBUFF3(*, 640)  - PE  INTEGER 

INBYT  - CU  INTEGER 
INDEX1  - CU  INTEGER 

INDEX 2 - CU  INTEGER 

INDEX3  - CU  INTEGER 

INDEX 4 - CU  INTEGER 

IPAGE  - CU  INTEGER 

LCH  - CU  LOGICAL 

LDIFFW  - CU  LOGICAL 

LF  - CU  LOGICAL 


- gives  the  byte  number  of  the  first 
byte  in  BUFFI. 

- controls  debug  print  out.  For 
routine  runs  equals  zero. 

- used  in  movement  of  data  within  time 
windows.  Difference  in  rows. 

- used  in  movement  of  data  within  time 
windows.  Difference  in  words. 

- frequency  currently  being  processed. 

- number  of  time  scans  in  input. 

- non-zero  if  a time  gap  was  found 
before  this  time  window. 

- user  supplied  factor  used  in  detecting 
glitches,  or  spikes,  in  the  data. 

- highest  frequency  processed. 

- input  buffer.  Equivalenced  to  NBUFF1. 

- output  buffer.  Equivalenced  to 
BUFF3(*,640) . 

- current  input  byte. 

- holds  intermediate  results  thru  out 
program. 

- holds  intermediate  results  thru  out 
program. 

- holds  intermediate  results  thru  out 
program. 

- holds  intermediate  results  thru  out 
program. 

- input  area  page  number  to  read  from 
next. 

- equivalenced  to  CH  to  facilitate 
shifting  and  masking. 

- equivalenced  to  DIFFW  to  facilitate 
shifting  and  masking. 

- equivalenced  to  F to  facilitiate 
shifting  and  masking. 
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LGAP  - CU  LOGICAL 


LINBYT  - CU  LOGICAL 


LNEW  - CU  LOGICAL 


LNGDCH  - CU  LOGICAL 


LNGDR  - CU  LOGICAL 


LNGDST  - CU  LOGICAL 


LNGT  - CU  LOGICAL 


LOFFSE  - CU  LOGICAL 


LOFREQ  - PE  INTEGER 
LOLD  - CU  LOGICAL 


LTSCAN  - CU  LOGICAL 


LT1  - CU  LOGICAL 


LT2  - CU  LOGICAL 


LT3  - CU  LOGICAL 


LT4  - CU  LOGICAL 


LT5  - CU  LOGICAL 


LT6-  CU  LOGICAL 


- equivalenced  to  GAP  to  facilitate 
shifting  and  masking. 

- equivalenced  to  INBYT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  NEW  to  facilitate 
shifting  and  masking. 

- equivalenced  to  NGDCH  to  facilitate 
shifting  and  masking. 

- equivalenced  to  NGDR  to  facilitate 
shifting  and  masking. 

- equivalenced  to  NGDST  to  facilitate 
shifting  and  masking. 

- equivalenced  to  NGT  to  facilitate 
shifting  and  masking. 

- equivalenced  to  OFFSET  to  facilitate 
shifting  and  masking. 

- lowest  frequency  processed. 

- equivalenced  to  OLD  to  facilitate 
shifting  and  masking. 

- equivalenced  to  TSCAN  to  facilitate 
shifting  and  masking. 

- equivalenced  to  T1  to  facilitate 
shifting  and  masking. 

- equivalenced  to  T2  to  facilitate 
shifting  and  masking. 

- equivalenced  to  T3  to  facilitate 
shifting  and  masking. 

- equivalenced  to  T4  to  facilitate 
shifting  and  masking. 

- equivalenced  to  T5  to  facilitate 
shifting  and  masking. 

- equivalenctd  to  T6  to  facilitate 
shifting  and  masking. 
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LT7  - CU  LOGICAL 

LTWSZR  - CU  LOGICAL 

NBUFF1(*,64)  - PE  INTEGER 

NCHAN  - CU  INTEGER 
NEW  - CU  INTEGER 

NGDCH  - CU  INTEGER 

NGDR  - CU  INTEGER 
NGDST  - CU  INTEGER 
NGT  - CU  INTEGER 
NROWS  - GU  INTEGER 

OFFSET  - CU  INTEGER 
OLD  - CU  INTEGER 

OPAGE  - CU  INTEGER 

OTIME(*)  - PE  INTEGER 

OVLAP  - CU  INTEGER 
PEN(*)  - PE  INTEGER 

PINT1(*)  - PE  INTEGER 
PINT2(*)  - PE  INTEGER 
PREAL1(*)  - PE  REAL 


- equivalenced  to  T7  to  facilitate 
shifting  and  masking. 

- equivalenced  to  TWSZu  to  facilitate 
shifting  and  masking. 

- input  buffer.  Also  referenced  by 

the  equivalenced  variable  RBUFF1 (4096) . 

- number  of  channels. 

- which  half  of  BUFF2  is  new.  Either 
one  or  two. 

- number  of  good  channels  after  variance 
check. 

- number  of  good  rows.  NGDCH  x TWSZR. 

- number  of  good  sites. 

- number  of  good  times.  NGDCH  x TWSZ. 

- number  of  rows  occupied  by  data  in 
BUFF2.  Equal  to  NSITE  x TWSZ/64. 

- used  in  index  calculation. 

- which  half  of  BUFF 2 is  old.  Either 
one  or  two. 

- page  in  output  area  to  be  written  to 
next. 

- time  in  deciseconds  from  the  beginning 
of  the  year  of  previous  data  word. 

- overlap  between  time  windows. 

- PE  number.  Initialized  in  block  area 
subroutine  to  be  1 in  pe  number  one, 

2 in  pe  number  two,  and  64  in  pe 
sixty  four. 

- used  for  intermediate  results  thru 
out  program. 

- used  for  intermediate  results  thru 
out  program. 

- used  for  intermediate  results  thru 
out  program. 


PREAL2 (*)  - PE  REAL 


- used  for  intermediate  results  thru 


RBUFF1(4096)  - PE  REAL 


RBUFF2(70400)  - PE  REAL 
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SITES (80)  - PE  INTEGER 

SITEGD(80)  - PE  INTEGER 

TIME(*)  - PE  INTEGER 

T0TSCN(*)  - PE  INTEGER 
TSCANS  - CU  INTEGER 

TVARFT(*)  - PE  REAL 

TWSZ  - CU  INTEGER 

TWSZR  - CU  INTEGER 

TWTIME(*)  - PE  INTEGER 

T1  - CU  INTEGER 
T2  - CU  INTEGER 
T3  - CU  INTEGER 
T4  - CU  INTEGER 


out  program. 

- input  buffer.  Equivalenced  to  NBUFFI 
and  IBUFFI. 

- intermediate  buffer  in  which  time 
windows  are  constructed.  Equivalenced 
to  BUFF2  and  ABUFF2. 

- give  physical  site  number  for  each 
site  used  for  output. 

- indicates  which  sites  passed  variance 
test . 

- time  in  deciseconds  from  the  beginning 
of  the  year  of  current  data  word. 

- number  of  time  scans  procesed  so  far. 

- number  of  time  scans  down  in  current 
time  windows. 

- holds  variance  factor  as  it  is 
modified  to  pass  50%  of  channels. 

- time  window  size.  Integer  power  of 
two  between  64  and  512. 

- number  of  rows  occupied  by  one  time 
window.  TSWZ/64, 

- time  in  deciseconds  from  the  beginning 
of  the  year  of  time  window  currently 
being  prepared  by  each  PE. 

- holds  intermediate  results  thru  out 
program. 

- holds  intermediate  results  thru  out 
program. 

- holds  intermediate  results  thru  out 
program. 

- holds  intermediate  results  thru  out 
program. 
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T5  - CU  INTEGER 


T6  - CU  INTEGER 


T7  - CU  INTEGER 


VARPT (*)  - PE  REAL 


- holds  intermediate  results  thru  out 
program. 

- holds  intermediate  results  thru  out 
program. 

- holds  intermediate  results  thru  out 
program 

- user  supplied  factor  used  in  detection 
of  dead  or  noisy  channels. 


LINKAGE 

In  addition  to  the  system  I/O  routines,  DEM2  requires  several  reloca- 
table modules  which  make  up  the  University  of  Illinois  FFT  routine.  These 
are: 


FFT.MAINREL 
FFT.SCRAMBLEREL 
FFT. TRAN REL 
FFT . CONSTANTSREL 


The  interface  to  these  modules  is  provided  by  the  subroutine  RUNFFT 
and  is  described  in  the  documentation  of  that  subroutine. 


DEM2  also  calls  subroutine  GTDATA,  C16T64,  C64T32,  and  ROWSUM  described 
in  the  section  on  subroutines. 


DISK  AREAS 


INDM2  - binary  input.  Described  in  DEMI  DISK  AREAS. 

0UTDM2  - binary  output.  Format: 

HEADER  - beginning  of  area. 

WORD  1:  Array  ID 

WORD2-1024 : zero 

TIME  WINDOW  - this  data  is  best  viewed  as  a two  dimensional  matrix. 

Each  column  will  be  wholly  contained  by  a PE  and  contains 
one  time  window.  There  are  64  columns,  one  per  PE. 

This  is  repeated  as  many  times  as  required  to  contain 
all  time  windows  processed.  The  format  of  each  column 
is : 
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WORD  1:  time  in  deciseconds  for  this  time  window. 

WORD  2:  number  of  good  channels. 

WORD  3-27:  physical  number  of  each  good  channel. 

WORD  28-640: 

First  Frequency  (channels  1-N) 


last  Frequency  (channels  1-N) 


CONPRM  - input  parameters  for  this  run,  in  binary. 

WORD1 : DEBUG.  Controls  debug  printout.  Normally  zero. 

WORD2:  time  window  size.  Integer  power  of  two  between  64  and  512. 

W0RD3:  over  lap.  Integer  between  zero  and  time  window  size 
minus  one. 

WORD4:  glitch  factor.  Normally  ten.  Floating  point. 

WORD5:  variance  factor.  Normally  ten.  Floating  point. 

WORD6:  low  frequency.  Integer  indicating  serial  number  of 

lowest  frequency. 

WORD7:  high  frequency.  Integer  indicating  serial  number  of 

highest  frequency. 

WORD8:  indicates  whether  vertical  or  horizontal  components  are 

to  be  processed.  1-vertical;  2-horizontal. 


D1SP2  - printed  output  consisting  of  a header  identifying  the  run,  informa- 
tion on  each  timing  error  located,  and  a trailer  indicating  normal 
termination . 
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DEM2  (tain  Prograa  Flowchart 


FKCOMB  ALGORITHM 


PURPOSE 

To  accomplish  rapid  signal  detection  using  frequency  beamforming,  i.e., 
frequency  wavenumber  analysis,  on  a large  number  of  beams. 

FUNCTIONAL  AND  THEORETICAL  DISCUSSION 

Frequency-wavenumber  (f-k)  spectral  estimation  is  a powerful  technique 
for  signal  detection  and  waveform  analysis  of  digitally  recorded  array  data. 
The  f-k  spectrum  of  a given  segment  of  array  output  is  the  squared  modulus 
of  the  multidimensional  Fourier  transform  of  the  data  in  time  and  space. 

The  f-k  spatial  representation  of  a propagating  wave  is  shown  schematically 
in  Figure  2.  Using  discrete  Fourier  analysis  in  the  frequency  domain,  the 
representation  can  be  thought  of  as  a series  of  layers  normal  to  the  frequency 
axis,  each  layer  representing  the  wavenumber  plane  at  a given  frequency.  The 
wave  is  thus  represented  as  a series  of  power  maxima  in  the  layers,  and  the 
locus  of  these  maxima  is  determined  by  the  phase  velocity  of  the  wave 
(Smart  1971) . 

The  advantage  of  this  process  is  that  propagating  wave  components  are 
easily  recognized  and  separated  from  one  another,  subject  to  the  limitations 
imposed  by  the  array  geometry,  sensor  weighting,  and  the  type  of  spectral 
smoothing  employed.  In  essence,  f-k  analysis  is  beamforming  in  the  fre- 
quency domain.  The  method  takes  advantage  of  the  fact  that  the  signal-to- 
noise  ratio  varies  with  frequency,  so  the  beamforming  is  done  frequency  by 
frequency.  Also  by  staying  in  the  frequency  domain  a great  many  beams  can 
be  examined  rapidly,  the  number  being  limited  only  by  the  resolution  cell  of 
the  array  response.  In  practice  this  means  that  the  azimuth  and  velocity 
of  a signal  need  not  be  assumed:  one  merely  accepts  the  beam  with  maximum 
power.  This  fact  is  important  for  signals  such  as  long-period  seismic 
surface  waves,  which  not  only  are  dispersive  (i.e.,  their  phase  velocities 
vary  with  frequency)  but  whose  arrival  azimuth  may  also  vary  with  frequency 
due  to  lateral  inhomo gene ties  in  their  paths. 
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Figure  2.  The  frequency  wavenumber  representation  of  a propagating 
wave . 


FKCOMB  is  a fast  f-k  analysis  program  that  was  used  in  an  automatic 
processing  system  for  microbarograph  array  data  (Smart  and  Flinn,  1971). 

Throughout  the  program  description  the  Fourier  transform  of  the  nth 

user-specified  channel  component  for  any  component  of  motion  at  angular 

iot 

frequency  w will  be  written  as  A^CuOe  n(co) . 

PROGRAM  DESCRIPTION 

FKCOMB  uses  f-k  analysis  for  continuous  processing  of  time-varying  data 
from  arrays  of  sensors.  Its  output  is  in  the  form  of  a bulletin  listing 
signal  detections  and  giving  for  each  the  phase  velocity,  back  azimuth, 
signal  power,  signal-to-noise  ratio,  and  F statistic  (related  to  signal- 
to-noise  ratio;  see  Smart  and  Flinn  1971)  as  a function  of  frequency  and 
arrival  time. 


FKCOMB  examines  time  windows  of  data  which  have  been  Fourier-transformed 
in  time  and  space.  Maxima  of  power  in  three-dimensional  f-k  space  are  auto- 
matically picked  and,  if  these  maxima  exceed  a specified  signal-to-noise 
detection  threshold  and  are  within  a specified  phase  velocity  range,  they 
are  listed  together  with  their  corresponding  back  azimuth,  phase  velocities, 
and  other  data.  (See  above.) 

There  are  also  two-dimensional  maxima,  which  are  places  that  are 
maximum  within  a given  wavenumber  plane  but  not  along  the  frequency  axis. 

If  such  two-dimensional  maxima  satisfy  the  specified  signal-to-noise  ratio 
and  phase  velocity  criteria,  and  if  their  corresponding  approximations  to 
group  velocity  (see  below)  are  reasonable,  these  maxima  are  also  listed  by 
the  processor. 

THE  FAST  f-k  ANALYSIS  ALGORITHM 

The  power  at  a given  frequency  and  wavenumber  is  computed  by  means  of 
equation  (1) : 


i\  1 1 v / c\  ia  (f)  2uik*r  i2 
P(f,k)  = b 1 A <f>e  n e nl 


'N  Vn 
n=l 


where 


f = frequency 
k = vector  wavenumber 

N = number  of  channel  components  for  the  component  of  motion 
n = channel  component  index 

r^  = vector  location  of  the  n’th  channel  component  with  respect 
to  an  arbitrary  origin 

An(f)ei0tn^  = Fourier  transform  of  the  n’th  channel  component 

A (f)  = amplitude  part  of  the  transform 
n 

eian(f)  _ phase  parJ;  Qf  the  transform  in  which  °n^  = the  phase  angle. 

Equation  (1)  is  evaluated  for  a matrix  of  wavenumber  values  at  a series 
of  discrete  frequencies,  as  specified  in  the  input  parameters;  it  can  be 
considered  as  a three-dimensional  space  with  frequency  being  one  dimension 
and  the  vector  wavenumber  k being  the  other  two  dimensions.  For  computa- 
tion, k is  resolved  into  a Cartesian  coordinate  system,  with  k related 

y 

to  geographic  north  and  k^  related  to  geographic  east. 

A wavenumber  value,  say  kQ,  is  related  to  the  phase  velocity  V by 
equation  (2)  : 

V = f/|kQ|. 

Stated  verbally,  phase  velocity  is  inversely  proportional  to  the  distance 
from  the  frequency  axis.  The  locus  of  constant  values  of  V is  a cone  in 
f-k  space,  with  the  apex  at  the  point  f = (k^k^)  = (0,0). 

For  f-k  analysis  the  power  is  calculated  in  a matrix  of  wavenumber 
values  separated  by  a grid  interval  Ak.  This  is  greatly  speeded  up  by 
using  the  relation  shown  in  equation  (3): 


A 'f)eian(f)e27Ti(k+k)*rn 


ia  (f)  277ik*r  2iTiAk*r 

= A (f)e  n e ne 

n 


n 
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Thus  if  a set  of  N terms  had  been  calculated  for  the  first  wavenumber  valuo 
k^,  the  values  at  = + Ak.  are  obtained  by  merely  multiplying  those 

terms  by  a factor  exp (+2TTiAk*  r ) . Therefore,  if  a regular  grid  is  used, 
only  one  set  of  kernels  exp(+2'.k_^*£n)  need  be  generated,  the  remaining 
values  being  obtained  with  successive  multiplication  by  the  invariant  ker- 
nels exp(+2niAk* r ). 

- — n 

THE  SEARCH  PROCESS 

The  frequency  wavenumber  search  using  triangular  and  rectangular  grids 
can  be  thought  of  as  taking  place  within  a cone  (Figure  3).  At  each 
frequency  searched,  the  process  takes  place  within  a search  disk  bounded 
by  the  intersection  of  the  cone  and  a constant  frequency  plane.  This 
search  disk  is  then  extended  by  a border  shaped  like  an  annulus  to  insure 
complete  coverage  by  the  grid.  Initially  a triangular  grid  is  used.  Once 
a maximum  is  found,  finer  square  grids  are  used,  utilizing  uphill  walks  from 
the  maximum  to  get  a better  estimate.  The  program  searches  from  lower  to 
upper  user-specified  frequencies. 

Beginning  at  the  point  of  greatest  power  on  the  coarse  grid,  the 
program  steps  out  in  a plane  of  constant  frequency  along  each  of  the  four 
cartesian  coordinate  directions  to  determine  the  direction  in  which  the 
power  is  rising,  and  it  continues  to  compute  successive  points  in  that 
direction  as  long  as  the  power  is  rising. 

When  the  power  begins  to  fall  off  in  the  direction  being  explored,  a 
new  direction  is  determined  and  followed,  and  the  process  is  repeated  until 
a place  is  reached  where  the  four  adjacent  points  in  f-k  space  all  show 
lower  power.  The  grid  spacing  is  then  reduced  by  a factor  of  6 and  the 
same  procedure  is  repeated  to  refine  the  location  of  the  power  maximum. 

The  amount  of  computation  required  is  about  an  order  of  magnitude  less 
than  would  be  required  for  computing  and  searching  the  complete  two- 
dimensional  spectral  matrix. 

All  two-dimensional  peaks  located  in  this  manner  are  then  checked  to 
see  if  they  are  also  maxima  in  the  frequency  direction  as  well;  such  peaks 
are  defined  by  the  equation: 
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Of  course  the  only  extrema  considered  are  the  maxima,  as  minima  are  of  no 
interest. 

The  program  checks  through  the  different  frequencies  as  follows: 
consider  the  power  in  the  wavenumber  plane  for  the  nth  frequency.  The  power 
maximum  in  this  layer  is  compared  with  the  power  and  location  of  the  maxi- 
mum in  the  (n-l)th  and  (n+l)th  layers;  for  clarity  of  exposition  we  call 


the  power  at  these  peaks  P , P , and  P 

n n-1  n+1 


If  both  P and  P , are 
n+1  n-1 


smaller  than  P^,  then  Pn  is  flagged  as  a three-dimensional  maximum.  If 
either  or  8reater  than  P^,  a check  is  made  on  the  respective 

positions  of  the  peaks  in  wavenumber  space;  suppose  that  P >P  . The 
vector  wavenumber  separation  between  the  location  of  P and  the  location 
Pn-1  calculated  and  the  following  logic  is  followed. 

(a)  If  the  separation  is  less  than  half  the  width  of  main  lobe  of 
the  array  response,  i.e.,  within  the  detection  cylinder,  the  two  peaks  are 
presumed  to  be  part  of  the  same  signal,  and  Pfi  is  not  designated  as  a three- 
dimensional  maximum. 

(b)  If  the  separation  is  greater  than  one  half  width  of  the  main  lobe 
of  the  array  response,  i.e.,  outside  of  the  detection  cylinder  P is 
presumed  to  be  part  of  another  signal.  To  see  where  P itself  is  nevertheless 
a three-dimensional  maximum,  the  power  is  checked  in  the  (n-l)th  and 

(n+1) th  layers  at  the  same  vector  wavenumber  as  the  peak  P : suppose  these 

are  and  respectively.  Then  Pn  > and  Pn  > Qn+1  This 

identification  process  is  continued  for  all  the  frequency  layers. 

Wavenumber  peaks  which  are  two-dimensional  but  not  three-dimensional, 
i.e.,  those  for  which 

3k  3k  u’  3f  r u 
x y 


may  nevertheless  appear  in  the  bulletin,  listed  separately  from  the  three- 
dimensional  peaks.  As  a check  on  the  validity  of  the  two-dimensional  peaks, 
a group  velocity  check  is  made  in  the  following  manner. 
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The  values  of  M/jAkj  are  calculated  for  each  peak  and  listed  in  the 
bulletin,  where  Af  is  the  frequency  difference  and  Ak  is  the  difference  in 
wavenumber  position  of  the  maxima  in  the  two  adjacent  layers.  The  ratio 
Af / | Ak | is  a first-order  approximation  to  the  group  velocity  of  a propagating 
signal;  thus,  if  a two-dimensional  peak  is  representative  of  the  power  in  a 
propagating  signal  at  that  given  frequency,  the  approximate  group  velocity 
should  have  a reasonable  value. 

A common  cause  of  spurious  two-dimensional  peaks  is  leakage  of  power 
(due  to  the  finite  length  of  the  time  window)  from  a large  peak  in  the 
spectrum.  It  has  previously  been  shown  (Smart,  1971)  that  this  leakage 
occurs  along  lines  of  constant  vavenumber  (i.e.,  Ak  = 0),  so  the  calculated 
approximation  to  group  velocity  Af / | Ak|  will  fall  outside  the  normal  range 
of  group  velocities  expected  for  seismic  signals  (2  to  3 kilometer/second). 

At  each  frequency,  only  the  first  large  two-dimensional  maximum 
is  considered  in  order  to  avoid  picking  up  power  entering  the  array  through 
sidelobes  in  the  array  response. 


J 


FKCOMB  at  present  picks  only  the  primary  peak  at  each  frequency. 
Experience  in  using  FKCOMB  to  analyze  long-period  seismic  data  has  shown 
that  nisnals  are  usually  detected  from  the  three-dimensional  maxima,  and  the 
main  interest  in  the  two-dimensional  peaks  is  to  provide  a more  complete 
spectrum  of  the  detected  signal  waveform. 


DATA  AREAS  AND  SYMBOL  DEFINITIONS 


By  default  in  CFD,  all  CU  variables  are  in  common  between  all  subrou- 
tines. All  PE  variables  have  been  put  in  the  common  block  FKMAIN.  All 
variables  may  thus  be  treated  as  a global.  This  is  the  principle  form  of 
communication  used  between  subroutines. 


ADJF(*)  - PE  INTEGER 


ADKX(4)  - PE  REAL 
ADKY(4)  - PE  REAL 


- used  to  hold  index  of  adjacent  fre- 
quency during  2-D  or  3-D  check. 

- used  in  computing  fine  grid. 

- used  in  computing  fine  grid. 
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ANGLE  - CU  REAL 


ARRAY  - CU  INTEGER 
AZ(*)  - PE  REAL 
BORDER  - CU  REAL 


CHANAV(*)  - PE  REAL 

CNTRL(*,6)  - PE  INTEGER 
COSDK(*)  - PE  REAL 
COUNT2 (*)  - PE  INTEGER 
COUNT3(*)  - PE  INTEGER 
DEBUG  - CU  INTEGER 
DELTAF  - CU  INTEGER 

DELTAY  - PE  REAL 
DELTAK  - PE  REAL 
DELTX (500)  - PE  REAL 

DELTY (500)  - PE  REAL 

DELX(*)  - PE  REAL 

DELY (*)  - PE  REAL 

DIST  - PE  REAL 


DKX  - CU  REAL 
DX(500)  - PE  REAL 

DY  (500)  - PE  REAL 


- input  parameter  indicating  angle  in 
degrees  to  be  searched.  Zero  means 
search  all  directions. 

- array  being  processed. 

- azimuth  of  signal  at  maximum  power. 

- width  of  border  of  search  disk  to 
insure  that  the  disk  is  adequately 
covered  by  the  coarse  grid. 

- used  in  calculation  of  signal  to 
noise  ratio. 

- constant  describing  each  array. 

- used  in  calculation  of  Re(Power) 

- number  of  2 dimensional  maximum  found. 

- number  of  3 dimensional  maximum  found. 

- controls  generation  of  debug  output. 

- difference  between  frequencies. 

SAMPLE  RATE /TIME  WINDOW  SIZE. 

- coarse  grid  spacing  in  Y direction. 

- coarse  grid  spacing  in  X direction. 

- X displacement  in  K-space  between 
successive  points  on  the  coarse  grid. 

- Y displacement  in  K-space  between 
successive  points  on  the  coarse  grid. 

- X-modification  in  K-space  for  direction 
being  searched  in  fine  grid  search. 

- Y-modif ication  in  K-space  for  direction 
being  searched  in  fine  grid  search. 

- distance  from  X-axis  (in  K-space)  to 
highest  point  on  the  coarse  grid  line 
being  processed. 

- user  input  rectangular  grid  spacing. 

- K-space  x-coordinate  of  each  point  ok 
coarse  grid. 

- K-space  y-coordinate  of  each  point  on 
coarse  grid. 
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FKX(* , 25)  - PE  REAL 


FFT(* ,612)  - PE  REAL 

FMAX  (*,25)  - PE  REAL 

FPMAX(*)  - PE  REAL 

FSTAT(*)  - PE  REAL 
HDKX  - CU  REAL 
HIFREQ  - CU  INTEGER 
I - CU  INTEGER 

IGO  - CU  INTEGER 

INBUF (* ,640)  - PE  INTEGER 

IND  - CU  INTEGER 

INDEX  - CU  INTEGER 

INFREQ  - CU  INTEGER 
IPOWER(* ,25)  - PE  REAL 

ITPOW(*)  - PE  REAL 
J - CU  INTEGER 
K(*)  - PE  REAL 
KERNEL (*,25)  - PE  REAL 

KSEP (*)  - PE  REAL 

KX  - PE  REAL 


- the  x-coordinate  in  K-space  of  the 
power  maximum  for  each  frequency, 

- equivalenced  to  INBUF91.28).  Holds 
FFT  output  from  DEM2. 

- the  maximum  power  for  each  freqeuncy 
processed. 

- power  maximum  for  frequency  currently 
being  processed. 

- fisher  statistic. 

- half  of  DKX. 

- highest  frequency  to  be  considered. 

- used  as  a "DO"  index  throughout 
FKCOMB. 

- used  to  calculate  point  arrangement  in 
lines  of  the  coarse  grid. 

- array  into  which  the  input  data  is 
read. 

- "DO"  index  in  subroutine  GRID  ranging 
from  second  to  last  coarse  grid  point. 

- indicates  frequency  layer  to  check 
when  determining  2 and  3-D  maxima. 

- frequency  currently  being  processed. 

- Im(Power)  of  each  channel  for  a point 
to  be  beamed  from. 

- Im(Power)  for  current  point. 

- used  as  a "DO"  index  throughout  FKCOMB. 

- used  in  calculation  of  group  velocity. 

- holds  intermediate  results  during 
power  calculation. 

- separation  in  K.-space  of  two  adjacent 
power  maxima. 

- X— coordinate  in  K-space  of  center  of 
search  disk. 


1 
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KY  - PE  REAL 

KXMAX(*)  - PE  REAL 

KYMAX(*)  - PE  REAL 

KXSEP (*)  - PE  REAL 

KYSEP(*)  - PE  REAL 

LINE  - CU  INTEGER 

LINEP1  - CU  INTEGER 

LINES  - CU  INTEGER 

LOCATE (*)  - PE  INTEGER 

LOC2D(*,25)  - PE  INTEGER 

LOC3D(8, 25)  - PE  INTEGER 

LOFREQ  - CU  INTEGER 
LOWER  - CU  REAL 

MNCHAN  - CU  INTEGER 
MODE 3 - CU  LOGICAL 

N - CU  INTEGER 

NGHAN(*)  - PE  INTEGER 

NFREQ  - CU  INTEGER 


tr  n in i n ~n n- iit niwilii iMtiliii 


'i —coordinate  in  K— space  of  center  of 
search  disk. 

- X— coordinate  in  K— space  of  current 
power  max. 

- Y— coordinate  in  K— space  oi  current 
power  max. 

- X part  of  separation  in  K-space  of  two 
adjacent  power  maxima. 

- Y part  of  separation  in  K-space  of  two 
adjacent  power  maxima. 

number  of  coarse  grid  lines  on  one  side 
of  the  vertical  axis. 

- equal  to  one  plus  the  number  of  lines 
on  the  coarse  grid. 

total  number  of  lines  on  coarse  grid 
equal  to  2*lines+l. 
location  of  point  on  the  coarse  grid 
having  maximum  power. 

• index  indicating  frequency  in  which 
each  max  was  found, 
index  indicating  frequency  in  which 
each  3-D  max  was  found. 

■ lowest  frequency  to  be  considered, 
user  input  limit  bound  on  phase 
velocity. 

maximum  number  of  channels  in  any  PE. 
mode  pattern  indicating  which  PE's 
have  detected  a 3-D  maximum, 
used  as  a "DO"  index  throughout 
FKCOMB. 

the  number  of  channels  of  data  for 
each  PE. 

number  of  frequencies  being  processed. 
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NMODE  - CU  INTEGER 

NPOINT  - CU  INTEGER 

NPTS( *)  - PE  INTEGER 
NTIMES  - CU  INTEGER 

OFFSET( *)  - PE  INTEGER 

PALE  - CU  INTEGER 
PERIOD (*)  - PE  REAL 
PINT1(*)  - PE  INTEGER 

PREAL1(*)  - PE  REAL 

PREAL2(*)  - PE  REAL 

RADIUS  - CU  REAL 
REFINE  - CU  INTEGER 

RINBUF(*,640)  - PE  REAL 
RPOWER(*,25)  - PE  REAL 

RTPOW( *)  - PE  REAL 
SAM  - CU  INTEGER 
SIGN  - CU  INTEGER 

SIGNAL (*)  - PE  REAL 
S INDK( *)  - PE  REAL 
SWITCH  - CU  INTEGER 

TPOWER(*)  - PE  REAL 


- mode  pattern  indicating  which  PE's  have 
ruled  out  a 3-D  maximum. 

- index  indicating  which  coarse  grid 
point  is  being  processed. 

- number  of  points  on  the  coarse  grid. 

- "DO"  index  controlling  refinement  loop 
in  subroutine  FNGRID. 

- used  in  address  calculation  for  access 
to  INBUF. 

- next  page  to  be  read  from  input  file. 

- period  of  signal  at  power  maximum. 

- holds  intermediate  results  throughout 
FKCOMB. 

- holds  intermediate  results  throughout 
FKCOMB. 

- holds  intermediate  results  throughout 
FKCOMB. 

- radius  of  search  disk. 

- input  parameter  indicating  number  of 
times  to  refine  fine  grid. 

- equivalenced  to  INBUF(1,1). 

- Re(Power)  of  each  channel  for  a point 
to  be  beamed  from. 

- Re(Power)  for  current  point. 

- sampling  rate. 

- indicates  sign  of  search  direction  on 
coarse  grid. 

- signal  to  noise  ratio. 

- used  in  calculation  of  Im(Power) . 

- passes  control  information  to  subrou- 
tine grid. 

- intermediate  power  maximum. 
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TWIN  - CU  INTEGER 
TWTIME(*)  - PE  INTEGER 


T1  - CU  INTEGER 

T2  - CU  INTEGER 

UPPER  - CU  REAL 

VEI. (*)  - PE  REAL 
X(*,25)  - PE  REAL 

XCOORD(*)  - PE  REAL 
Y (*,25)  - PE  REAL 

YCOORD(*)  - PE  REAL 
YMAX (50)  - PE  REAL 

YILMI  - CU  INTEGER 


YPOINT(50)  - PE  REAL 
YTOP  - CU  INTEGER 


- time  window  length. 

- time  in  deciseconds  from  beginning  of 
year  of  beginning  of  time  window  being 
processed . 

- scratch  variable  used  throughout 
FKCOMB. 

- temporary  variable  used  throughout 
FKCOMB. 

- user  input  limit  bound  on  phase 
velocity. 

- velocity  of  signal  at  maximum  power. 

- X-coordinate  of  channels  used  for 
each  PE. 

- X-coordinates  of  all  channels. 

- Y-coordinate  of  channels  used  for  each 
PE. 

- Y-coordinates  of  all  channels. 

- maximum  Y displacement  for  each 
coarse  grid  line. 

- one  less  than  number  of  coarse  grid 
points  on  coarse  grid  line  under 
construction. 

- number  of  points  on  each  coarse  grid 
line. 

- maximum  vertical  offset  of  coarse  grid 
line  being  processed. 


LINKAGE 

FKCOMB  utilizes  the  standard  COS,  SIN,  and  SORT  functions  supplied 
with  CFD.  It  also  utilizes  the  I/O  package  supplied  by  Institute  for 
Advanced  Computation  (IAC). 


FKCOMB  also  calls  subroutine  GRID,  FNGRID,  CHECKR,  OUTPUT,  REALE,  IMG, 
and  MAX  which  are  discussed  in  the  section  on  subroutine. 
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DISK  AREAS 


CONPRM 
ST CORD 
FKIN 
FKDISP 


Contains  user  supplied  run  time  parameters. 

Contains  coordinates  of  sensors  of  array  being  processed. 
Main  input  area  for  FKCOMB.  See  output  of  DEM2. 

Area  for  printed  output  from  FKCOMB. 
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The  preliminary  version  of  FKCOMB  demonstrates  the  feasibility  of 
processing  seismic  data  on  ILLIAC.  The  features  and  restrictions  of  this 
version  were  planned  so  that  the  code  could  be  brought  up  in  a minimum 
amount  of  time  and  still  both  demonstrate  the  feasibility  of  putting  this 
and  other  algorithms  on  ILLIAC  and  provide  a solid  base  for  development  of 
a version  for  routine  processing  of  long  period  data.  The  code  has  been 
debugged  to  test  sample  data.  Processing  long  period  data  routinely  would 
require  further  debugging  of  the  code  and  removal  of  the  following  restric- 
tions : 

/Arrays  - The  current  version  has  been  tested  on  LASA  data.  Coding  is 
included  to  handle  ALPA  and  NORSAR  data,  but  has  not  been  tested. 
Processing  these  arrays  differs  only  in  D0il  and  DEM2;  the  FKCOMB 
analysis  is  identical. 

Overlap  - Overlapping  of  time  windows  is  coded  but  untested  in  the  current 
version.  This  capability  requires  the  routing  of  data  from  one  time 
window  to  the  next  and  was  not  tested  in  the  preliminary  version. 

Motion  Components  - The  Gene  Smart  version  of  FKCOMB  allows  for  processing 
of  horizontal  as  well  as  vertical  motion  components.  In  actual  use, 
vertical  components  are  of  primary  interest.  The  processing  of 
horizontal  components  was  therefore  designed  into  the  algorithm,  but 
never  coded.  If  found  to  be  necessary,  this  option  can  be  added 
without  major  revision. 

Channels  - The  ILLIAC  version  currently  allows  the  user  to  specify  the  first 
channel  and  the  last  channel  in  the  tape  description  in  the  block  data 
subroutine  initialization  of  the  array  CNTRL.  This  was  thought  to  be 
sufficient.  In  practice  it  was  found  that  it  is  often  necessary  to 
ignore  arbitrary  channels.  This  addition  can  be  added  quite  simply  by 
zeroing  out  any  channel  that  is  not  used  so  that  it  does  not  affect  the 
power  maximum.  It  is  still  necessary  to  keep  track  of  the  number  of 
real  channels  in  t"der  to  correctly  calculate  average  power  and  other 


Variance  (detection  of  dead  and  noisy  channels)  - As  originally  designed, 
after  detection,  dead  and  noisy  channels  were  eliminated  by  actually 
removing  them  from  the  data  stream.  During  implementation,  especially 
after  consideration  of  the  removal  of  unwanted  channels  (see  above), 
it  was  decided  that  a much  better  action  after  the  detection  of  a dead 
or  noisy  channel  is  to  zero  it  out  and  treat  it  exactly  as  an  unwanted 
channel.  Thus,  the  code  for  removal  of  the  channel  was  never  debugged 
and  currently,  though  dead  and  noisy  channels  are  detected,  they  are 
not  removed.  This  addition  requires  the  same  effort  as  removing 
unwanted  channels. 

Stripping  - Stripping  is  not  included  in  the  current  version  of  FKCOMB,  but 
is  a straight  forward  addition.  It  is  an  option  in  Lhe  Smart  version. 

Half  Frequencies  - The  Gene  Smart  version  of  FKCOMB  zero  fills  before  FFT 
in  order  to  obtain  twice  as  many  frequencies.  The  added  freqeuncy 
resolution  gained  is  minimal,  so  this  feature  was  not  included  in  the 
ILL1AC  version.  In  addition  to  the  loss  of  some  frequency  resolution, 
there  are  several  possible  side  affects  which  may  result  from  this 
difference  in  algorithm  between  the  two  versions.  In  particular,  the 
loss  of  resolution  is  most  important  for  short  time  window  lengths. 
Discrepancies  between  the  results  of  the  two  versions  must  be  analyzed 
before  determining  whether  this  feature  should  be  added  to  the  ILLIAC 
version.  If  found  necessary,  the  addition  is  straight  forward. 

Timing  - The  run  time  is  reducible  by  at  least  an  order  of  magnitude  by 
local  optimization  of  code.  Much  of  this  consists  of  the  insertion 
of  assembly  language  code  at  critical  points.  The  repeated  calculation 
of  exponentials  in  the  analysis  portion  of  the  routine  can  be  removed 
with  little  work  and  would  result  in  a significant  increase  in  speed. 
Due  to  these  considerations,  the  lack  of  anything  but  wall  clock  time 
from  ILLIAC,  and  the  small  amount  of  data  currently  being  processed, 
it  it  impossible  to  give  any  timing  data  other  than  the  theoretical 
time  given  elsewhere  in  this  report. 
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The  ILLIAC  version  of  FKCOMB  was  tested  by  the  comparison  of  it  with 
the  Gene  Smart  version  currently  available  on  the  UCLA  IBM  360/91.  A 
seventeen  minute  and  four  second  tune  period  starting  at  1628  GMT  plus 
12.0  seconds,  day  79,  1972  was  chosen  for  testing  due  to  the  existence  of 
a large  signal  in  the  short  period  data  from  that  time  period.  The  following 
modifications  were  made  to  a standard  deck  used  to  run  FKCOMB  on  the  UCLA 
360/91: 

1.  SLOP  was  changed  from  4.0  to  1000.0  so  that  dead  and  noisy  channels 
would  not  be  eliminated  so  as  to  duplicate  the  conditions  of  the  ILLIAC 
version. 

2.  Time  window  size  and  INC  were  both  set  to  256  to  provide  for  a 256- 
point  time  window  and  no  overlap. 

3.  Channel  numbers  and  coordinates  were  modified  so  as  to  include 
channels  5 and  6,  which  normally  are  not  used.  This  was  necessary  due  to 

the  restriction  in  the  ILLIAC  version  discussed  earlier.  Channel  coordinates 
of  zero  were  used  for  these  channels,  as  they  are  not  valid  long  period 
channels. 

The  seismograms  of  an  event  in  Honshu,  Japan  (t=15  57  50.4,  40. 8N, 

141. 9E,  March  19,  1972,  ^=6)  chosen  to  compare  the  ILLIAC  version  of  FKCOMB 
to  the  version  programmed  for  the  IBM  360/91  are  shown  in  Figure  4.  The 
position  of  four  nonoverlapping  256  second  time  windows  analyzed  are  also 
indicated  on  the  figure.  The  time  windows  cover  the  initial  portion  of  the 
Rayleigh  wave  train  also  including  the  arrival  times  of  some  seismic  core 
phases  (PKKP  and  PKKS).  The  results  of  parallel  runs  of  the  two  versions  are 
given  in  Figure  5.  In  the  figure  the  results  of  the  ILLIAC  version  are  listed 
in  the  columns  marked  14  and  the  results  of  the  Smart  version  on  the  360/91 
are  listed  in  the  columns  marked  91.  The  dashes  indicate  some  events  located 
by  the  ILLIAC  version  of  FKCOMB  were  not  listed  by  the  Staart  version  of  FKCOMB. 
The  Smart  version  eliminates  two  dimensional  maxima  which  do  not  lie  along 
physically  reasonable  phase  velocity  lines  through  f-k  space.  This  avoids 
reporting  erroneous  detections  which  are  artifacts  of  the  f-k  spectral 
process.  The  sorting  algorithm  for  review  of  results  from  FKCOMB  analysis 
on  ILLIAC  will  also  eliminate  these  detections  but  was  not  implemented  in  the 


preliminary  version.  The  result  agree  well  in  general,  the  differences  can  be 
attributed  mainly  to  differences  in  the  details  of  search  routines  used.  The 
ve3  ^cities  associated  with  the  maximum  F statistics  in  the  tables  are  reasonably 
close  to  the  phase  velocities  of  Rayleigh  waves  at  the  corresponding  periods, 
and  the  azimuths  to  those  of  the  direction  of  the  epicenter.  Some  of  the 
values  seem  to  be  somewhat  high  but  it  must  be  kept  in  mind  that  the  relatively 
small  size  of  the  LASA  array  does  not  permit  the  precise  determination  of 
phase  velocities;  thus,  inaccuracies  in  phase  velocities  can  be  expected.  The 
results  indicate  that  the  two  programmed  versions  of  FKCOMB  work  properly  and 
yield  comparable  results. 
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Figure  5.  Sample  FK.C0MB  Output 
Sheet  2 of  2 


C16T64 


DESCRIPTION 

L16I64  converts  16  bit  raw  data  as  sent  from  the  seismic  arrays  to  64  bit 
II.LIAC  floating  point  format.  The  format  of  the  16  bit  data  depends  upon  the 
value  of  ARRAY.  If  ARRA5  is  equal  to  1 or  2,  the  data  consists  of  a 14  hit 
two's  complement  integer  left  justified  in  a 16-bit  byte.  The  data  is 
extracted  as  shown: 


M 


13,14  15 


VALUE  = S x(-213)+M 


If  ARRAY  is  equal  to  3,  the  data  consists  of  a 12  bit  two's  complement 
number  and  a 4 bit  gain  code,  all  right  justified.  The  data  is  extracted 
as  follows: 


d 

M 

0 

3 4 

15 

VALUE  = (S  x(-2U)+M)  x 210_G 


Most  of  the  code  consist  of  inserted  assembly  language  code  in  a CFD 
driver.  If  ARRAY  is  not  1,  2 or  3,  an  error  message  is  generated. 


ARGUMENTS 


1 N ( * ) - PE  REAL  - input  data. 
OUT ( *)  - PE  REAL  - output  data. 
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C16T64 


CUT  3 2 


DESCRIPTION 

C64I32  creates  a 32-bit  ILLIAC  floating  point  number  from  a 64-bit 
floating  point  number.  The  32-bit  result  is  left  in  the  outer  half  of  a 64 
bit  word.  The  sign  is  transferred  directly.  The  mantissa  is  truncated 

without  rounding  to  24  bits.  The  new  exponent  is  computed  via  the  following 
formula : 

NE  = new  exponent 
OE  = old  exponent 
NE  = (OE-16384)+64 

No  check  is  made  to  see  if  the  input  number  is  outside  the  range  of 
32-bit  representation.  The  inner  half  of  the  word  is  zero. 

ARGUMENTS 

IN(*)  - PE  REAL  - 64-bit  floating  point  input. 

OUT(*)  - PE  REAL  - 32-bit  floating  point  output. 
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DESCRIPTION 

Subroutine  CHECKR  examines  the  power  maximums  found  at  each  frequency 
to  determine  which  maximums  are  three  dimensional  and  which  are  two  dimen- 
sional in  frequency  wavenumber  space. 

ARGUMENTS 

IPMAX(*)  - PE  REAL  - power  maximum  for  frequency  being  processed. 
LOCATE(*)  - PE  INTEGER  - location  of  point. 

LOC2D(*,25)  - PE  INTEGER  - index  indicating  frequency  in  which  each  3-D 

maximum  was  found. 

LOC3D(*,25)  - PE  INTEGER  - index  indicating  frequency  in  which  each  2-D 

maximum  was  found. 


j 


■ 

' 

I 
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FPI1AX(I)>  FPMAX ( I + 1 ) 


CNVTIM 


DESCRIPTION 

CNVTIM  is  called  when  the  next  item  to  be  retrieved  is  the  timing  data. 
The  input  pointer  is  pointing  at  the  first  byte  of  the  timing  data.  The 
format  of  the  timing  data  is  dependent  upon  the  site.  Word  2 of  CNTRL(*, 
ARRAY)  give  the  format.  Currently  there  are  two  formats  accepted.  If  CNTRL 
(2, ARRAY)  is  equal  to  1,  the  time  is  in  the  form  of  deciseconds  from  the 
beginning  of  the  year  and  takes  up  2 bytes.  No  conversion  is  necessary. 

The  byte  is  retrieved  and  stored  in  TIME.  If  CNTRL(2 .ARRAY)  equals  2,  the 
timing  data  is  in  4 bit  characters  giving  DDDHHMMSS.  Three  16-bit  bytes  are 
retrieved  and  TIME  calculated  to  be  in  deciseconds  from  the  beginning  of  the 
year.  If  CNTRL (2. ARRAY)  has  a value  other  than  1 or  2.  an  error  message 
is  generated. 

ARGUMENTS 

CNTRL (2, ARRAY)  - PE  INTEGER  - value  indicates  format  of  timing  data. 
TIME(*)  - PE  INTEGER  - time  in  deciseconds  from  the  beginning  of  the 

year  as  calculated  by  CNVTIM. 
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DESCRIPTION 

FNGRID  refines  the  value  of  die  power  maximum  found  on  the  coarse  grid. 
It  computes  power  in  the  east,  south  west  and  north  directions  about  the  point 
found  on  the  coarse  grid  until  the  power  does  not  increase  in  any  direction. 

A square  grid  is  imposed  around  the  point  with  a spacing  between  points 
which  is  one  sixth  (1/6)  the  spacing  on  the  coarse  grid.  Figure  6 shoys  the 
square  grid  about  p.  the  max  on  the  coarse  grid.  The  grid  is  refined  the 
number  of  times  specified  by  the  user. 

ARGUMENTS 

REFINE  - CU  INTEGER  - number  of  times  to  refine  grid. 

DKX  - CU  REAL  - rectangular  grid  spacing. 

KXMAX(*)  - PE  REAL  - X-coordinate  in  K.  space  of  current  power  maximum. 

KYMAX(*)  - PE  REAL  - Y-coordinate  in  K space  of  current  power  maximum. 

FPMAX ( *)  - PE  REAL  - power  maximum  for  frequency  being  processed. 


-57- 


FNGRID 


GETBYT 


DESCRIPTION 

GETBYT  retrieves  the  16  bit  byte  pointed  at  by  INBYT  from  the  disk  area 
INPUT.  Two  levels  of  buffering  are  used,  ADBBUF  in  ADB  and  INBUF  in  core. 
First,  a check  is  made  to  see  if  the  address  wanted  is  currently  in  ADB.  If 
so,  the  byte  is  retrieved.  Otherwise,  the  ADB  buffer  must  be  refilled  from 
core.  A similar  check  is  made  to  see  if  the  address  is  in  core.  If  so, 

ADB  is  refilled  and  the  byte  retrieved.  Otherwise,  core  is  refilled  from  the 
appropriate  location  on  disk,  ADB  refilled,  and  the  byte  retrieved. 
Initialization  is  such  that  on  the  first  call,  all  buffers  are  marked  empty 
and  will  be  filled  before  the  byte  is  retrieved.  As  long  as  bytes  requested 
are  consecutive,  the  refill  of  buffers  will  be  minimal.  Most  bytes  actually 
retrieved  by  DEMI  are  consecutive,  with  the  extra  logic  necessary  to  handle 
skips  between  time  scans.  Accesses  are  assumed  to  be  always  increasing,  that 
is  the  address  requested  is  assumed  to  be  greater  than  the  previous  one. 

ARGUMENTS 

INPTB  - CU  INTEGER  - serial  number  or  address  within  disk  area  of  byte 

desired. 

INBYT  - CU  INTEGER  - returns  byte  address  by  INPTB. 
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DESCRIPTION 

Grid  computes  the  radius  of  the  search  disk  for  every  frequency  and  the 
spacing  between  grid  points  in  the  X and  Y direction.  It  passes  the  KX  and 
KY  coordinate  for  the  first  point  on  the  grid  and  the  change  in  X and  Y 
directions  to  move  progressively  along  the  grid.  On  the  second  call  to  grid 
the  actual  KX  and  KY  coordinates  for  the  maximum  power  location  found  on  the 
coarse  grid  are  returned. 

As  mentioned  above,  a coarse  triangular  wavenumber  grid  is  chosen 
initially  in  order  to  keep  computing  time  to  a minimum.  However,  the 
spacing  must  be  fine  enough  to  ensure  that  at  least  one  point  falls  on  a 
main  lobe  above  the  height  of  the  largest  side  lobe  of  the  array  response. 
For  example,  if  the  largest  side  lobe  of  the  array  response  has  a peak  9 dB 
down  from  the  maximum  of  the  response,  the  grid  spacing  has  to  be  selected 
so  that  at  least  one  point  falls  wit » in  the  9 dB  contour  of  the  main  lobe 
for  any  signal  in  the  velocity  range  of  interest. 

In  the  initial  coarse-grid  computation  of  the  wavenumber  spectrum, 
FKCOMB  employs  an  equilateral  triangular  grid  because  it  requires  fewer 
points  and  less  computation  time  than  does  a square  grid.  Nevertheless, 
the  convention  has  been  adopted  that  users  will  specify  the  optimum  square 
grid  interval  required  for  the  given  sensor  array  and  FKCOMB  will  convert 
that  to  the  corresponding  equilateral  triangular  grid  spacing  (Figure  7). 
figure  S shows  the  point  arrangement  on  the  coarse  grid. 

ARGUMENTS 

BORDER  - CU  REAL  - width  of  border  of  search  disk. 

DKX  - CU  REAL  - rectangular  grid  spacing. 

ANGLE  - CU  REAL  - radius  of  search  disk. 

LOWER  - CU  REAL  - lower  limit  of  phase  velocity. 

UPPER  - CU  REAL  - upper  limit  of  phase  velocity 

DELTX(500)  - PE  REAL  - grid  spacing  on  Y axis 

DELTY(500)  - PE  REAL  - grid  spacing  on  X axis 

KXMAX(-)  - PE  REAL  - X coordinate  in  K space  of  current  power  maximum. 

KYMAX(*)  - PE  REAL  - Y coordinate  in  K space  of  current  power  maximum. 
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UKX  is  the  grid-spacing  input  parameter  calculated  from  the  array  response 
AX  is  the  spacing  between  points  on  the  coarse  grid  in  X direction 
AY  is  the  spacing  between  points  on  the  coarse  grid  in  the  Y direction. 

Figure  7.  Coarse  Grid  Spacing. 
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Figure  8.  Point  Arrangement  on  Coarse  Grid. 
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GTDATA 


DESCRIPTION 

GTDATA  retrieves  bytes  sequentially  from  the  disk  area  INDM2  (IN  DEM2) . 
The  input  data  goes  through  two  levels  of  buffering,  one  in  ADB,  one  in  core. 
Each  call  on  GTDATA  results  in  a 16-bit  byte  being  placed,  right  :ustified, 
in  INBYT,  and  the  updating  of  pointers  to  the  next  16-bit  byte.  The  pointers 
are  initialized  in  the  main  driver,  DEM 2 , to  reflect  the  fact  that  all  buffers 
are  empty  before  the  first  call. 

ARGUMENTS 

INBYT  - CU  INTEGER  - returns  the  next  16-bit  byte  right  justified  and 

padded  with  zeroes. 





GT  DAT  A 


^ BYTE  ^ 
IN  CURRENT 
N\f)B  WORD‘S 


X IS  X 

BYTE  IN  ADB 
\ BUFFER 


^ IS  \ 

BYTE  IN 
LORE 


REFILL  CORE 
FRO'I  DISK 


REFILL  ADB 
BUFFER  FROil 
CORE 


RETURN 
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DESCRIPTION 

IMG  takes  a complex  number  composed  of  a real  and  an  imaginary  part 
stored  in  the  outer  and  inner  32-bit  portions  of  a 64-bit  wood,  extracts  the 
imaginary  part,  which  is  the  inner  part  in  ILLIAC  terminology,  and  converts 
it  to  ILLIAC  64-bit  floating  point  format.  The  routine  consists  almost  en- 
tirely of  assembly  language  code  inserted  in  a CFD  driver. 

The  complex  number  is  in  the  following  format: 
outer  sign  inner  sign 


Outer  Mantissa 


34  40  63 


7 8 9 


Inner  Mantissa 


15  16 


Bit  numbers 


The  result  is  in  the  f>  lowing  format: 


sign  bit 


15  16 


The  sign  bit  is  transferred  without  alteration.  The  mantissa  is  expanded 
from  24-bits  to  48-bits  by  padding  24-bits  of  zero  on  the  right.  The 
exponent  is  calculated  via  the  following  formula: 

NE  = new  exponent. 

E = old  exponent. 

NE  = (E-64)+16384. 

ARGUMENTS 


IN(*)  - PE  REAL  - input  complex  number. 

0UT(*)  - PE  REAL  - output  64-bit  floating  point. 
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DESCRIPTION 

MAX  finds  the  maximum  of  64  fixed  point  numbers  one  per  PE,  and  returns 
that  value  in  each  PE.  The  method  used  is  to  compare  each  PE's  value  to  that 
1,  2,  4,  8,  16,  and  32  away  via  the  route  instruction  and  to  substitute  the 
greater  value  each  time.  The  mode  is  assumed  to  be  ON  upon  entry.  The  inner 
branch  is  accomplished  by  temporarily  inhibiting  selective  PE's.  Subroutine 
MAX  is  used  to  determine  the  maximum  value  of  a variable  in  each  PE  when 
initializing  loops  in  the  code.  For  example  if  a loop  control  is  needed  for 
the  number  of  channels  MAX  determines  the  maximum  number  of  channels  stored 
in  the  PE's.  The  MODE  is  disabled  in  a given  PE  when  the  number  of  channels 
stored  in  it  is  reached.  The  MODE  will  be  disabled  in  all  PE's  when  the 
loop  control  equals  the  maximum  number  of  channels  found  by  subroutine  MAX. 

ARGUMENTS 

I(*)  - PE  INTEGER  - input  parameter  whose  maximum  is  desired. 

MAX(*)  - PE  INTEGER  - maximum  value  of  !(*).  Identical  in  all  PE's. 
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OUTPUT 


DESCRI PTLON 

The  statistics  associated  with  each  power  maximum  are  computed  in  sub- 
routine output.  The  fisher  statistic,  signal  to  noise  ratio,  azimuth, 
velocity  and  period  are  output  for  each  maximum  all  stored  for  each  the  time 
window.  See  flowchart  for  output. 

ARGUMENTS 

FPMAX(*)  - PE  REAL  - power  maximum  for  freojency  being  processed. 

KXMAX(,,:)  - PE  REAL  - X-coordinate  in  K space  of  current  power  maximum. 

KYMAX ( *)  - PE  REAL  - Y-coordinate  in  K space  of  current  power  maximum. 

AZ(*)  - PE  REAL  - velocity  of  signal  at  maximum  power. 

VEL(*)  - PE  REAL  - azimuth  of  signal  at  maximum  power. 

FSTAT ( *)  - PE  REAL  - fisher  statistic  at  maximum  power. 

SIGNAL (*)  - PE  REAL  - signal-to-noise  ratio  at  maximum  power. 

PERIOD( *)  - PE  REAL  - period  of  signal  at  power  maximum. 
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PUTBYT 


DESCRIPTION 

PUTBYT  isolates  the  right  hand  16  bits  of  OUBYT  and  inserts  them  in  the 
output  file  indicated  by  ARRAY.  The  byte  goes  through  two  levels  of  buffering, 
a one  word  ADB  buffer,  ADBOUT  and  a 4096  word  buffer  (one  for  each  output  file) 
OUTBUF.  Since  there  is  only  one  ADB  buffer,  it  is  necessary  to  save  and 
restore  it  whenever  output  arrays  are  changed.  It  is  also  necessary  to  empty 
the  buffers  when  processing  is  completed.  These  functions  are  carried  out 
by  the  main  driver,  DEMI. 

It  is  required  that  the  mbst  recent  output  always  be  available  in  core. 

For  this  reason,  PUTBYT  never  empties  the  core  buffers,  but  writes  the  first 
75%  to  disk  and  moves  the  last  25%  to  the  beginning  of  the  buffer  and  posi- 
tions the  pointers  appropriately. 


ARGUMENTS 


OUBYT  - CU  INTEGER  - holds  the  byte  to  be  output. 
ARRAY  - CU  INTEGER  - array  being  processed. 
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RDPRM 


DESCRIPTION 

RDPRM  is  provided  to  input  by  parameters  necessary  for  DEMI  processing. 
Currently,  the  only  parameter  necessary  is  DEBUG,  which  controls  debug 
printout.  Since  this  is  the  only  parameter,  it  is  currently  compiled  in  and 
no  input  from  disk  is  necessary.  Future  expansion  may  require  the  addition  of 
code  to  facilitate  reading  parameters  from  disk. 

ARGUMENTS 


DEBUG  - CU  INTEGER  - set  to  zero  for  routine  processing. 


REALE 


DESCRIPTION 

REA,,E  takes  a compLex  number  composed  of  a real  and  an  imaginary  part 
stored  in  the  inner  and  outer  32-bit  portions  of  a 64-bit  word,  extracts  the 
real  part,  wh.ch  is  the  outer  part  in  ILLIAC  terminology,  and  converts  it  to 
iLLIAC  64— bit  floating  point  format.  The  routine  consists  almost  completely 
of  assembly  language  code  inserted  in  a CFD  driver. 


The  complex  number  is  in  the  following  format: 

outer  inner 

sign  sign 


bit  numbers 


The  result  is  in  the  following  format: 

sign 


The  sign  bit  is  transferred  without  alteration.  The  mantissa  is  expanded 
from  24  bits  to  48  bits  by  padding  24  bi's  of  zero  on  the  right.  The  exponent 
is  calculated  via  the  following  formula: 

NE  = new  exponent. 

E = old  exponent. 

NE  = (E-64) +16384 


ARGUMENTS 

INf*)  - PE  REAL  - input  complex  number. 

OUT(*)  - PE  REAL  - output  64-bit  floating  point. 
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ROWSUM 


DESCRIPTION 

Function  ROWSUM  returns  the  sum  of  64  floating  point  numbers,  one  per  PE. 
The  result  is  identical  in  all  PE's  except  for  differences  caused  by  adding 
the  numbers  in  different  orders.  The  algorithm  used  is  the  "standard  ILL  4C 
logsum"  method  whereby  partial  sums  are  compiled  for  each  sequence  of  2,  4,  8, 
16,  32,  and  then  all  64  PE's  by  routing  by  1,  2,  4,  8,  16,  and  32.  The 
routine  is  coded  via  insered  ASK  statements  within  a CFD  driver.  It  is 
assumed  that  MODE=ON  when  ROWSUM  is  called. 

ARGUMENTS 


R(*)  - PE  REAL  - the  numbers  to  be  summed. 

ROWSUM(*)  - PE  REAL  - the  result,  replicated  in  all  PE's. 
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RUNFFT 


DESCRIPTION 

RUNIFT  is  a CFD  driver  which  makes  the  University  of  Illinois  FFT 
routine  accessible  to  a CFD  program  as  a normal  subroutine  call.  It  is 
primarily  assembly  language  code  to  facilitate  the  difference  in  parameter 
linkage.  The  FFT  routine  is  called  with  the  following  parameters: 

time  window  size  - TOSZ 
number  of  time  windows  - NGDCH 
input  data  - BUFF2 (* , 1 ,OLD) 

The  FFT  routine  requires  that  the  parameters  be  pointed  to  by  acar  2. 
Acar  2 contains  the  address  of  a 3 word  block.  Each  word  of  this  block 
contains  the  address  of  one  of  the  arguments.  Word  2 give  the  number  of 
discrete  time  windows.  The  third  word  contains  the  address  of  the  data  to 
be  FFT'ed.  The  University  of  Illinois  routine  currently  leaves  the  machine 
in  32-bit  mode.  The  driver  returns  the  machine  to  64-bit  mode  before  pro- 
cessing continues. 

ARC 'MEN  TS 

The  arguments  to  FFT  are  passed  ir,  common  via  AD8 . The  variables 
involved  are: 

TWSZ  - CU  INTEGER  - power  of  2 between  64  and  312  specifying  time 

window  size. 

NGDCH  - CU  INTEGER  - number  of  time  windows,  to  FFT. 

BUFF2(*,1,0LD)  - PE  REAL  - oasses  input  data  and  output  data. 


.'%■{»  jrf" 
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DATA  PREPARATION 

Two  input  files  must  be  supplied  by  the  user  at  run  time.  The  first, 

RUN. INPUT,  contains  the  raw  input  as  read  from  a long  period  tape  This  data 
]s  currently  provided  by  reading  the  tape  on  the  360/44  at  SDAC  and  shipping 
the  resulting  data  file  to  the  ILLIAC  site  over  the  ARPANET'  via  FTP. 

The  second  file,  RUN.PARAM,  contains  run-time  parameters  in  binary  form. 
There  is  currently  no  straight  forward  means  for  creating  such  a file.  During 
testing  of  the  program,  the  file  was  created  by  assembling  and  running  a job 
on  ILLIAC  which  writes  the  data  to  a disk  file.  The  parameters  are  assembled 
in  so  that  the  assembler  converts  the  data  from  symbolic  form  to  binary  form. 
The  source  of  the  assembly  language  program  currently  used  is: 

BEGIN 
START : : 

PUT  CON'PRM , CONPRMMEM , f OMPLWD : 

PAUSE  CCMPLWD; 

HALT; 

HALT; 

CONPRM : AREA  "CONPPM" ; 

COMPLWD: WDS  1; 

BLK; 

CONPREMMEM: DATA 

0,  % DEBUG=0 

64,  % TIME  WINDOW  LENTTH=64 

0,  Z OVERLAPS 

10-0,  % GLITCH  FACTOR=10.G 

10.0  % VARIANCE  FACTOR=10.0 

2»  % FREQUENCY  2 IS  LOWEST  FREQUENCY 

9>  % FREQUENCY  9 IS  HIGHEST  FREQUENCY 

0,  % 0 MEANS  VERTICAL  MOTION 

0,  % THIS  WORD  IS  NOT  USED 

0,  % THIS  WORD  IS  NOT  USED 


0, 

% THIS  WORD  IS  NOT  USED 

. 0J4 , 

% DKX 

2.90, 

% LOWER  VELOCITY  LIMIT. 

9999999. , 

1 UPPER  VELOCITY  LIMIT 

0.0, 

% ANGLE=0 .0 ; SEARCH  IN  ALL 

DIRECTIONS. 

2, 

1 REFINE  FINE  GRID  TWICE 

1. 

».  SAMPLING  RATE  IS  ONE  PER 

SECOND 

0; 

END  START. 

/ END  OF  PARAMETERS. 

The  primary  input  file  used  to  run  this  job,  assuming  thut  the  above 
source  is  contained  in  the  file  PA RAM. ASK  is: 

ASK  PARAM. ASK, PARAM. REL , PARAM. LIST 
LINKED 

ISV  PARAM. ISV 
INCLUDE  PARAM. REL 
END 

ALLOC  FKCOMB .MAPFILE,  1 
RUN  PARAM. ISV 

MOVE  14DM:C0NPRM,  RUN. PARAM 
DALLOC  1 


PROGRAM  EXECUTION 

In  order  to  run  the  ILLIAC  Long-period  Seismic  Array  Processing  System 
the  following  files  must  be  stored  on  the  ILLIAe  central  file  system: 

DEMI.  ISV 
DEM2.  ISV 
FKCOMB. ISV 
FKCOMB.  MAPFILE 
RUN. INPUT 
RUN. COORDINATES 
RUN. PIF 

RUN . PARAMETERS . 
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All  uf  tlie  above  files  except  KUN.  PARAMETERS  and  RUN. INPUT  are  provided. 
RUN. PARAMETERS  contains  user  specified  run  time  parameters  and  should  be 
prepared  tor  each  run.  RUN. INPUT  contains  the  long-period  seismic  data.  See 
Data  preparation,  section  Iv.b  for  description  of  these  files. 

RUN. PI F (primary  input  file)  contains  the  job  control  language  tiiat 
follows  : 

RUN. PIP 

ALLOC  FKCOMB . MAPI?  I L E , 1 
MOVE  RUN. COORDINATES , 14DM:STC0RD 
MOVE  R U N . P A RAM , 1 4 DM : CO N P RM 
MOVE  RUN.  INPUT, 14DM: INPUT 
RUN  DEMI . ISV, MAXTIM=600 
RUN  DEM2 . ISV,MAXT1M=600 
RUN  FKCOMB. 1SV,MAXT1M=600 
MOVE  14  DM : FKD I S P , RUN . OUTPUT 
DALLOC  1. 

The  command  ACL  "SUBMIT  RUN.PIF"  enters  the  job  onto  the  ILLIAC  batch 
queue.  For  details  of  job  submittal  and  monitoring  see  ACL  USER's  GUIDE. 

Upon  completion  the  disk  Files  RUN. OUTPUT  and  RUN.POF  will  have  been  created. 
RUN.POF  contains  the  system  responses  to  the  job  command.  RUN. OUTPUT  is  an 
8-bit  ASCII,  file  suitable  for  printing  containing  bulletins  for  all  2-D  and 
3-D  maximums  found  and  the  associated  statistics. 
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ACAR  - any  one  of  four  accumulators  in  the  CU. 

AUB  - Advast  Data  Buffer.  64  words  of  fast  access  memory  used  by  the  CU. 

ASK  - refers  to  the  ILMAC  assembler  and  the  assembly  language. 

CFD  - a Fortran  like  language  developed  by  NASA  which  compiles  on  an  IBM  360 
and  outputs  ILL1AC  assembly  language. 

Components  of  motion  - the  vertical,  north,  or  east  direction  of  wave  motion. 

Cl  - control  unit.  Hie  ILLIAC  Control  Unit  decodes  the  instruction  stream 

to  be  executed  by  the  PE's.  It  also  executes  a limited  set  of  instructions. 

Data  record  - a number  od  data  scans.  If  not  every  data  scan  has  a time  word, 
the  data  record  will  contain  a time  word  corresponding  to  the  time  of 
the  first  scan. 

Data  scan  - a unit  containing  the  values  for  all  the  data  channels  at  a seismic 
station  for  a given  instant,  which  may  also  be  on  the  scan.  During  data 
retrieval  one  scan  at  a time  is  processed. 

MODE  - a pattern  of  64-bits  which  controls  the  write  protection  of  registers 
and  memory  for  each  PE. 

PE  - process.  *<*  element.  One  of  64  ILLIAC  processors  which  execute  the  same 
instructs  stream  in  lock  step. 

Route  instruction  - provides  an  efficient  means  for  communication  between  PE's. 

Time  word  - word  holding  the  time  of  the  first  data  scan  in  a u^a  frame. 

Because  the  data  scans  are  spaced  equally  in  time,  the  time  ol  any  scan 
in  the  frame  is  known  if  the  time  of  the  first  scan  is  known.  GMT  is 
used  in  all  time  words. 
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